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Abstract 

Packet-dispersion based measurement tools insert pairs of probe 
packets with a known separation into the network for transmission 
over a unicast path or a multicast tree. Samples of the separation 
between the probe pairs at the destination(s) are observed. Heuristic 
techniques are then used by these tools to estimate the path charac- 
teristics from the observations. In this paper we present a queueing 
theoretic setting for packet-dispersion based probing. Analogous to 
network tomography, we develop techniques to estimate the parame- 
ters of the arrival process to the individual links from the samples of 
the output separations, i.e., from the end-to-end measurements. 

The links are modeled as independent discrete time queues with 
i.i.d. arrivals. We first obtain an algorithm to obtain the (joint) dis- 
tribution of the separation between the probes at the destination(s) 
for a given distribution of the spacing at the input. The parameter 
estimates of the arrival process are obtained as the minimizer of a cost 
function between the empirical and calculated distributions. We also 
carry out extensive simulations and numerical experiments to study 
the performance of the estimation algorithm under the fairly 'harsh' 
conditions of non stationarity of the arrival process. We find that the 
estimations work fairly well for two queues in series and for multicast. 

Keywords: Network tomography, packet-pair measurements, discrete 
time queues, packet-pair tomography, multicast queues. 
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1 Introduction and Background 



Network tomography is the estimation of detailed statistics of performance 
parameters from aggregate or end-to-end measurements of some measurable 
quantity. The term was first coined in pQ, where the problem of estimat- 
ing flow volumes from link volumes was analyzed. This is also called the 
traffic matrix estimation problem and has been extensively studied in the 
transportation literature, for which, an excellent survey is available in [2J. 
The traffic estimation problem is now being addressed by the networking 
community ( [31 SI El El E] • Another network tomography problem is that of 
estimating link delay statistics from the samples of end-to-end delays expe- 
rienced by multicast packets [8, 9j. Much of the work in this area has been 
in estimating the link or segment delay distributions from the end-to-end 
measurements. See [10] for a fairly comprehensive survey of both delay and 
traffic estimation tomography problems and solutions related to networking. 

The bandwidth on a network path has also been of significant interest 
and a number of tools have been developed for its estimation. Two kinds 
of bandwidth related metrics are often estimated. The bottleneck bandwidth 
is the maximum transfer rate that could be achieved over a path and is 
essentially the minimum of the transmission rates on the links in the path. 
The available bandwidth is the portion of the bottleneck bandwidth not 
used by competing traffic and depends on the traffic load at the inputs to 
the links. 

Many tools have been developed to measure the bottleneck and available 
bandwidths. To the best of our knowledge, all these tools are based on either 
'packet-pair' or 'packet-dispersion' techniques. The packet-pair technique 
was first described in [llj to measure the bottleneck bandwidth on a path 
where the links are rate-servers. In this technique two back-to-back packets 
of equal length are transmitted by the source. It can be shown that the ratio 
of the probe packet length to the separation between them at the receiver 
is the service rate for the packet at the bottleneck link on the path. This 
idea has been used to develop many tools, e.g., pathchar [12] and clink 
|13j . that measure the bottleneck capacity of a network path. Inspired by 
the packet-pair technique, the packet-dispersion, also called packet-spacing, 
technique has been developed to measure the available bandwidth of an 
Internet path [HI [HI [161 El US]- I n this technique, a source transmits 
a number of probe packets with a predetermined separation, say d. The 
samples of the separation at the receiver, say d r , are then used to estimate 
the available bandwidth on the path, pathload |14j and pathChirp |15j are 
some examples of tools that use the packet-spacing technique to measure 
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the available bandwidth. See [16] for an excellent survey of the different 
packet-spacing techniques used in the measurement tools and [T7] for an 
experimental comparison of the tools. A non cooperative version that does 
not require a measurement-enabled receiver is described in [18j . 

Packet-pair and packet-dispersion techniques have an important advan- 
tage in that they do not require the sender and the receiver to be synchro- 
nized in time. This makes them important in practice. However, the tools 
that employ these techniques are based on informal arguments. There have 
been excellent statistical studies based on experiments on packet-pair tech- 
niques (e.g., [19, 20J). Recently, there have been a few theoretical mod- 
els developed for packet dispersion based bandwidth probing techniques 
[3H |221 E31 |231 Eg ESI EZ|. In [2H [22] it is shown that the output dis- 
persion of probe pair samples the sum of three correlated random processes 
that are derived from the traffic arrival process. The asymptotics of these 
and the probing process are also analysed. The measurement methods of 
some bandwidth estimation tools are then related to the analytical models 
developed. These results are then extended to multihop paths in |23j . In 
[24] . the network is assumed to be a time-invariant, min-plus system that 
has an unknown service curve. This service curve estimation method is 
developed for different probing techniques. 

While the above works provide interesting theoretical insights into the 
probing process, our work is more closely related to that of |27} \25\ 126] . In 
[27] , a probe train is assumed and the actual delay of each of the probe 
packets is assumed to be known at the output of the path. The delay 
sequence is assumed to form a Markov chain whose transition probability 
matrix is estimated. An inversion is defined to estimate the characteristics 
of the arrival process from the transition probabilities of the delay Markov 
chain. In [25], the bottleneck link is modeled as an M/D/l queue, and the 
expected dispersion is calculated using the transient analysis. The method 
of moments is then used to estimate the arrival rate. A similar method 
is followed in [26J where the more general M/G/l queue is considered and 
Takacs' integro-differential equation is used to obtain the distribution of the 
dispersion. 

In this paper we develop techniques to estimate the parameters of the 
arrival process distribution on the individual links on a path or a multicast 
tree from packet-dispersion samples in a queuing theoretic setting. Packet- 
pair probes with predetermined separation between them are injected at a 
source. These probe packets traverse discrete-time queues with i.i.d. packet 
arrivals on a path or in a multicast tree. Samples of the separation between 
the probe packets at the output nodes of the path or the multicast tree 
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are used to estimate the parameters of the arrival process distribution at 
the individual links. The approach is as follows. First, for a single queue 
and for a given separation between the probes at the input, we derive the 
conditional distribution of the separation between the probes at the output 
of the queue in terms of the distribution of the arrival process to the discrete 
time queue. The key intermediate result here is the joint distribution of 
the number arrivals to the queue and the number of departures from the 
queue between the slots in which the probes are injected. Then the output 
separation is obtained for any given distribution of the input separation. 
Since the input separation distribution is known, (we assume it to be fixed) 
this is applied recursively on the path to obtain the separation distribution 
at the outputs. The parameter estimator is the minimizer of a suitable 
distance function between the empirical and the theoretical distributions of 
the output separations. 

We first develop our results by assuming that the probe packets are 
served after all the packets that arrive in their slot, i.e., the probe packets 
have the least priority among the packets that arrive in their slot. Extension 
to the case of the probe packet having the same priority (or any arbitrary 
fixed priority) as the other packets that arrive in the slot will be along the 
same lines and is also presented. We present numerical results only for the 
first case, i.e., when the probes have lowest priority. 

Our modeling assumption is similar to that of [27] in that we also con- 
sider discrete time queues. While [27] considers the delay process of the 
probe sequence, we consider the dispersion of the packet pair. The inversion 
methods are also different. While we consider a general i.i.d. arrival pro- 
cess for the discrete time queeu, [25] considers an M/D/l queue and obtains 
approximations, [26J considers an M/G/l queue. 

The rest of the paper is organized as follows. In the next section, we 
describe the problem setting and the notation that we will be using in the 
paper. In Section[3]we consider a single pair of probes that are injected into 
a stationary discrete-time queue and develop an algorithm for computing 
the distribution of the separation between a pair of probes at the output 
of a single queue when the distribution of the separation at the input of 
the queue and also the arrival process distribution are known. In Section HI 
we consider sending the probes over a multi-link path and over a multicast 
tree and obtain the distribution of the output separation(s) between the 
probes at their outputs. In Section [5] we describe the method to estimate 
the parameters of the arrival process at the queues. The numerical evalua- 
tion of the estimates for the case of Poisson arrivals are discussed in detail 
in Section [H In Section [7J we describe a generalization and an extension 
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Figure 1: The convention for A n , X n and Q 



to the basic theory developed in Section In particular, we outline the 
method for computing the output separation distribution when the probes 
have equal priority with the packets arriving in the same slot. We conclude 
with a discussion in Section [8] where we explore connections with an ear- 
lier theoretical study of packet-dispersion techniques. We also contrast our 
results with delay-tomography results. 



2 Notation and Preliminaries 

We first describe the notation used in this paper. For pedagogical purposes, 
we consider a packet queuing system. In this section we will consider a 
single-server queue with infinite buffers and FCFS service discipline. The 
service time of each packet is equal to the slot length. The convention 
will be as follows. Packet arrivals to the queue in slot n will occur at the 
beginning of the slot and the departures will occur at the end of the slot. 
Thus arrivals in a slot will be available for departure in the slot. We observe 
the queue at the beginning of a slot, just after the arrival instant. Let 
A n denote the number of packets that arrive in slot n, X n the number of 
departures in slot n and Q n the number of packets in the queue when the 
queue is observed at the beginning of the slot. This convention is shown in 
Figure [TJ Of course, X n can only be or 1. Further, for to < n, A m ^ n : — 
A m +A m+ i + - ■ -+A n will denote the cumulative arrivals in the interval [to, n] 
and X m ^ n := X m + X m+ i + • • • + X n will denote the cumulative departures 
in the interval [m, n]. In particular, X m ^ m = X m and A m ^ m = A m . Note that 
X m ^ n is a function of Q m and [A m , . . . , A n ]. Xm,n will denote the number 
of departures in slots [in, n] that the same arrival sequence [A m ,...,A n ] 
would cause with Q m = q. Observe that this makes I^n independent of the 
distribution of Q m and a function of only [A m , . . . , A n \. The work-conserving, 
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FCFS service discipline will be applied at the granularity of the batch of 
packets that arrive in a slot and the packets in a batch will be queued in a 
random order. The queue evolution equation as per the above convention 
will be 



Qn+l — Qn + A n +i X n , 
Qn+m = Qn ~\~ -^-n+l,n+m -^n,n+ra— 1 



In the rest of this paper we assume that the arrival sequence A n is 
i.i.d. and let pk '■= Pr (^4o = k) . Recall that if A := E (Aq) < 1 then the 
queue length process is stationary and the moment generating function Q(z) 
of the stationary queue length distribution 7r q := Pr (Qq = q) is given by 



(l-\)(z-l) V { z ) 
z - V{z) 



where V(z) is the moment generating function of p^- In fact 7r q can be 
obtained using the following recursion. 



7T« 



Pr(Qo = 9) 
l-E(A)) = 1-A 
^(to(I-Po)) 



^ (7T,_l(l -Pi) -TT q - 2 P2 ~ 
~ Kip q -l - TlQPq-l) 



for q = 0, 
for q = 1, 

for q > 1. 



(2) 



3 Passage of a Probe Pair through a Discrete Time 
Queue 

In this section we consider injecting a single probe pair into a stationary 
queue and derive the distribution of the output separation, d , as a function 
of the input separation, dj, and the parameters of the packet arrival process 
{^4 n }. For a given distribution of the input separation, the distribution of 
d a is obtained as 

oo 

Pr (d = s) = Pr {do = s\di = j) Pr (d t = j) . (3) 
i=i 

Let us consider a queue in steady state. Consider two probe packets, 
denoted by P\ and P2, inserted into the discrete time queue in slots and 
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di respectively. Since we are considering stationary queues, without loss of 
generality, we can relabel the slots in this manner. These probe packets are 
enqueued along with the other packets that arrive in their respective slots 
and are available for transmission in the same slot in which they arrive. In 
this section, we assume that the probe packets have the least priority among 
the packets that arrive in their slot. 

Let D\ and D 2 denote the departure slots for the first and the second 
probe respectively. Let d Q := D 2 — D± denote the dispersion of the probe 
packets at the output of the queue. Clearly, 

D l = Q , 

D 2 = di + Q dl = di + Q + l + A 1)di - X 0) di~i , 
d :=D 2 -D 1 = di + l + A 1A - X 0A _i . (4) 

d Q — di is the packet- dispersion. In the above, Qo does not include probe 
packet P\ and Qd t is the number of packets in the queue at the beginning of 
slot di not including P 2 . Also, A di does not include P 2 while X^^-i includes 
the possible departure of Pi. Note that Xq^-i > -^0 = 1 since the queue 
has Qo + 1 packets (including Pi) in slot 0. Observe that d Q , the separation 
between the probes at the output of the queue, is affected by the arrivals 
and departures between the arrivals of Pi and P 2 , i.e., by A\^ and Xq^-i- 
Thus, to obtain the conditional distribution of the output separation, given 
the input separation, we need to obtain the joint distribution of the number 
of arrivals and departures between the two probe pairs. Computation of 
this joint distribution is the key to the computation of the distribution of 
the output separation. 

Prom Eqn. [H for a given di, say di = m, the distribution of the output 
separation d Q can be expressed as 

m 

Pr (d = s\di = m) = Pr (X , m -i = /, A 1>m = s - m + I - 1) (5) 

1=1 

for s = 1,2, Since we are assuming that we are probing a stationary 

queue, the first probe packet Pi will see the queue in steady state. We can 
then obtain the joint distribution of Xo )m _i and Ai^ m as 

Pr (^o,m-i = A\ yrn = j) 

00 

= ^ Pr (X , m -i = I, A 1)Tn = j\Q = q) Pr (Q = q) 

00 

= Pr ( X oVl = *> A hm = j) (Qo = q) ■ (6) 

q=0 



7 



We can see that Qq and A\^ m are independent of each other, but Xo jm _i 
is dependent on both Qq and on the finite sequence A\,... ,A m . We make 
the reasonable assumption that the distribution of Qq is known since it can 
be derived from the distribution of Aq, albeit as a function of its param- 
eters. Hence, to obtain the distribution of d , we first need to obtain the 
joint distribution of Xo jm _i and A\^ m conditioned on Qq = q, i.e., the joint 
distribution of Xq m _ 1 and A^ m . This is derived in the following subsection. 



3.1 Joint Distribution of A im and X? 



0, rn- 



In this section we first obtain a recursion on m, for the joint distribution 
of A\^ m and Xq m _ 1 and then, for arbitrary q, express the joint distribution 
of A\^ m and X^ m _ 1 in terms of the joint distribution of Ai >m and X^ m _ l . 
We thus begin with the case of Qq = 0. From our notation, this means that 
when the first probe packet is injected into the queue in slot it is the only 
packet in the queue. Since there is the probe packet to be transmitted in 
slot 0, X® = 1. Thus we have the following base equations for the joint 
distribution. 



Pr (X% = l,A ltl =j) 




j) for 1 = 1, 
otherwise. 



(7) 



Now, for arbitrary m > 1, we can write the following one-step recursion 
which also spells out the kind of recursion we seek. 

= J2 J2 Pr ( X tn-1 = h Ai, m = j\X^ m _ 2 = I, Ai, m _i = n) 

i n 

x Pr (X m _ 2 = i, A ljm _! = n) . (8) 

For ease of exposition, define 

4>] := Pr {Xl^ = I, A hm = j\Xl m _ 2 = i, A hm ^ = n) . 

We now consider three cases. 

1. If XQ m _ 2 = i < I — 1 or ^o,m-2 = i > I, then = 0. This is because, 
for there to be / departures in slots [0, m — 1], the number of departures 
in slots 0, . . . , (m — 2), denoted by i, should be either I or I — 1. 
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2. If lj m _ 2 = i = I — 1, then there has to be a departure in slot m — 1. 
Since Qo = 0, this is possible if, and only if, A\^ m _i = n > I — 1 
because only then will the total number of arrivals in slots [0, m — 1], 
including Pi, be greater than or equal to I. Hence, for n < I — 1, 
ifil'™ = and for Ai :in _i = n > I — 1, lj m _[ = / is a sure event and 
Al t m = j requires that j — n packets have to arrive in slot m giving us 

= Pr (An = j 

3. For lj m _ 2 = i = I, there should be no departure in slot in — 1. 
Since Qo = 0, P± has departed in slot 0. Hence the total number 
of departures in slots [1, m — 2], i.e., / — 1, should be the same as 
Ai tTn _i = n. Therefore, for n 7^ I — 1, i/y - = 0. Thus for i = I, only 

n = I — 1 will yield a non zero ibl'"" when the number of arrivals in 
slot m is j — n = j — I + 1, i.e., for i = / and n = / — 1, we have 
^ = Pr(A m =j-l + l). 

Combining the above three cases, we can rewrite the recursion in Eqn.[5]for 
the joint distribution of X$ m _ 1 and Ai^ m as 

Pr (Xl m _ x = l,A 1>m =j) 
3 

= Pr ( A m = J - n) Pr (X ° m -2 = I ~ 1, Msn-x = n) 
n=l-l 

+ Pr(A m =j-l + 1) Pr (X° m _ 2 = Z, Ai jtn _i = i - l) . (9) 

We obtain the joint distribution Xq j and Ai^ m by a suitable trans- 
formation of the joint distribution of Xg m _ 1 and A\ m derived above. To 

obtain Pr \ X^ m _ 1 = l,Ai jTn = j\ we consider the following cases. Clearly, 

this probability is zero if I > m (the number of departures is greater than 
the number of slots) and for j < 0. 

1. The number of departures in slots [0, m — 1] is less than the number 
of slots, i.e., I < m. We consider two sub cases. 

(a) Clearly, q > I — 1 is not possible because in that case, the number 
of departures including the probe packet Pi would have been at 
least Z + l. Therefore, Pr (x£ m _i = l,A 1:m = j\ = for q>l-l. 

(b) For q < I — 1, the event 

{(Ai, . . . , A m ) I XQ m _ 1 = 1} 
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is the same as the event 

{(A 1 ,...,A m )\X° 0tm _ 1 =l-q} 

and so 

{(Ai,...,4n) | Xl m _ 1 = l,A 1<m =j} = 
{(A 1 ,...,A m ) \ Xj m _j = / - g,Ai, m = j} 

Thus, for ^ < I — 1, we get 

Pr(^m-i = J 'V = i) = 

Pr (Xj ro _ 1 = i-g,Ai )ro =j) 

2. If I = m, then there is a departure in every slot in [0, m — 1]. We once 
again consider two sub cases. 

(a) If Qo = q < {I — 1), then the event 

{(^i,...,A m ) | X« m _ 1 = m,A ltm =j} 
is the same as the union of the disjoint events 

{(A u ...,A m ) | Xj ro _! = i, Ai, m = j} 
for f = m — q, . . . , m. Hence, for this case, we can write 

Pr (Km-l = I.V = j) = 

£ Pr(X ° im _ 1 = t,A 1 , m = i). 

t=m— q 

(b) If the queue starts with Qo = 9 > — 1) packets, then there will 
be a departure in every slot and X^ m _ 1 = I is a sure event and 

so for this case Pr (x^ m _ 1 = /, Ai )Tn = j) = Pr (A 1>m = j). 
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Summarizing the above, we have the following transformation. 



'Pr(X° >m _ 1 =l-q,A 1}m =j) 

for I < m and q < I — 1 

for / = m and g < I — 1 

Pr (Ai >m = j) 

for / = m and g > Z — 1 
otherwise. 



(10) 



Once the joint distribution of -Xgm-l an< ^ * s found using the joint 
distribution of X$ m _ 1 and Ai^ m , the joint distribution of Xo jm _i and Ai^ m , 
can be found from Eqn. [6l 

The infinite summations in Eqns. [6] and [3] will need to be truncated in 
practice. Let N\ and N2 respectively denote the truncation limits. Also, 
in Eqn. [5] we need to limit the range of the output separation for which 
we will compute the probabilities. Let ./V3 be this limit. The algorithm for 
computing the distribution of the output separation is summarized in Fig. [2l 
The recursion in step 3b of Fig. [2] is computationally the most intensive part 
of the algorithm. It can be shown that this step requires 

( N 3 (N 3 + 1) 5^ A^ 2 (jV 2 + 1) _ N 2 (N 2 + l)(2iV 2 + 1) 



V" 
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additions and multiplications. 



4 Passage of a Probe Pair through a Network of 
Independent Queues 

We now characterize the output separation between the probes as they pass through 
multiple queues. Since we are able to express the distribution of the output separa- 
tion of the probe pair in terms of the distributions of the input separation and the 
arrival process distribution, we can obtain the distribution of the output separation 
at the output of a path of independent queues fairly easily. 

We first consider the case when the probe pair passes through a path consisting 
of K independent queues as shown in Figure [3l The probe packets leaving queue k 
are placed in queue (fc + 1) after an arbitrary but fixed delay. Since the delays 
are fixed, the separation between the probe pairs at the input to queue (fc + 1) is 
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1. Input: Arrival distribution {pk}k>o 

Input separation distribution {Pr (di = fc)}fc>i. 



2. Obtain n q using the following equations. 

!1 - E (Aq) = 1 - A for q = 0, 

^(tto(I-Po)) for 9 = 1, 

^ (7T 9 _l(l - pi) - 7T 9 _2]32 7TlPg_l - 7T p 9 -l) for Q > 1. 

3. Obtain the joint distribution of {Xq m _ 1 , Ai. m } for each q using the following 
equations. 



(a) Base Equations 

Pr (Xl = l,A 1A =j) 

(b) Recursions 



Pr (At = j) for I = 1 and j > 
otherwise 



Pr (*o° m -i = J, A hm = j) = £ Pr (A m = j - n) Pr (X ° ro _ 2 = 1-1, Ax. m -\ = n) 



n=l-l 



(c) Transformations 



PlW,m-l =! ' A l.m=j) = 



+ Pr (A m = j - / + 1) Pr (X ° ro _ 2 = Z, = I - l) . 



'Pr (X$ m _ 1 = I - q, A\. m = j) I < m and q < I — 1 

Pr (*0,m-l = *i ^i,m = j) Z = m and g < / - 1 

Pr (Ai iTO = j) Z — m and q > I — 1 

otherwise 



4. Obtain the joint distribution of {Xo im _i, Aim} by using the following equa- 
tions. 

Pr (X , ro _l = Z, A ljm = j) = J2 Pr (^0,m-l = i. M,m = j) *T, 

5. Obtain the output separation distribution for a given input separation, di, 
say di = to. 

m 

Pr (cZ G = s\d, = m) = ^ Pr (X , m _i = Z, Ai >m = s - to + Z - 1) for s = 1, 2, . . . , N 3 . 



6. Obtain the distribution of the outjput separation for a given input separation 
distribution. 



N 2 



Pr (d„ = s)=J2 Pr ( d o = s \ d i = j) p r {<k = J) ■ 



Ficmrp 9: Alcrnrithm t.n enmnutp thp Histrihntinn nf tbp nutnnt. SFvnara.tinn 
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do 



di 



d K 



Queue 1 



Queue 2 



Queue 3 



Figure 3: K independent queues and the path of the probe sequentially 
through the K queues. 



d, 



Queue 1 



Queue 2 



j 



Queue 3 



Figure 4: A three-queue multicast tree. 



the same as that at the output of queue k. (Generalization to variable delay is 
presented in Section [7]) Let do be the separation at the input to queue 1 and dk 
the separation between the probe pairs at the output of queue k with d out :— dx- 
For k = 1, 2, . . . , K, let A„ denote the i.i.d. arrival sequence into queue k, p[ k \ 

A lk) (k) 

the distribution of Aq and 7r] the stationary distribution of the queue occupancy 
at the beginning of a slot. Probe packets Pi and P2 are injected into queue 1 in 
slots and do respectively with a known distribution of do. The distribution of dk 
is obtained from the recursive relation 

00 

Pr (d k = s) = Pr (4 = s|4-i = 3) Pr (dfc-i - j) ■ (11) 

3=1 

The distribution of dx can be obtained by repeated use of the algorithm in the 
previous section from the distribution of do- Notice that, if do = 1? then d± has the 
same distribution as A^\ 

The second case that is of interest is when the probe pairs traverse a multicast 
tree of independent queues. Figure 2] shows an example of a three-node tree. Here 
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copies of the probe packet are made at the output of Queue 1 and the copies are 
simultaneously placed in queues 2 and 3 after arbitrary but fixed delays on each 
of the links. As with the case of a path, since the delays are fixed, the separation 
between the probe pairs at the input to queues 2 and 3 is the same as that at the 
output of queue 1. We are interested in the joint distribution of the separation at 
the output of both queues 2 and 3, d 2 and d 3 respectively. This joint distribution 
is given by 

Pr (d 2 — r, d 3 — s) 

oc 

= Pr ( rf 2 = r,d 3 = s\d! = t) Pr (dj = t) 
t=i 

oo 

= ^ Pr ( d 2 = r\dt = t) Pr (d 3 = s\d! = t) Pr (rfi = t) (12) 

t=i 

The last equality follows because conditional on d\, d 2 and d 3 are independent. The 
individual conditional probabilities can be obtained as before. 

The above ideas can be extended to compute the joint distribution of the sep- 
aration at the leaves of an arbitrary tree rooted at the node which introduces the 
probe pairs into the network. 

5 Estimating the Arrival Process Parameters 

The algorithm developed for computing the distribution of the output separation 
between the probes can, in principle, be used to estimate the traffic parameters in 
the network provided the parameters are identifiable from the distribution. How- 
ever, often the parameters may not be identifiable from the distribution, i.e., the 
same distribution may be resulted from two different sets of traffic parameters. For 
instanct, consider a series of two queues with independent Poisson traffic with means 
Ai and A2. The output separation has the same distribution for (Ai,A2) = (0, A) 
and (Ai,A2) = (A, 0). We performed some simulation experiments to judge the 
viability of using this distribution for traffic parameter estimation. However, these 
simulations use ideal settings like independent Poisson traffic arrivals to different 
queues which are usually not strictly satisfied by practical networks, and as a result, 
should be taken only as proof of concept. They do not guarantee that the method 
will work in practical networks. 

Consider a rooted tree of K independent queues. Probe pairs are introduced 
into this network at the root of the tree with unit separation. While the theory 
allows for an arbitrary but known distribution of the input separation, we choose a 
fixed input separation of one slot. This minimizes the computational requirements. 

The separation between the probe pairs at the output of the leaves is observed. 
We assume that there is no loss of probe pairs in the network. We also assume that 
the link delays, i.e., the delay between the output of a queue and the input to the 
next queue on a path is arbitrary but fixed. 
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(k) 

For k = 1, . . . K , let A„ be the arrival sequence of packets into queue k. For 
queue k, is an i.i.d. sequence and further, the sequences for the different queues 

(k) 

are independent. Let A& be the vector of parameters of the distribution of A and 
let us define A := [Ai, . . . , Xk]- Our goal is to estimate A by observing only the 
separation of the probes at the output of the leaves of the paths. Let ji, j 2 , ■ . ■ ,j K 
denote these leaf nodes and dj 1 , dj 2 , . . . , dj K the separation of the probe pairs at the 
output of these queues. 

We send a sequence of probe pairs and corresponding to every probe pair we 
obtain an observation of the vector d := [dj 1 , dj 2 , . . . , dj K ]. From these observations 
we can obtain the empirical distribution of d, denoted by 'I'd- We probe with a very 
low rate, i.e., the interval between the probe pairs is made very large so that the 
probing load is negligible. This allows us to assume that the stationary distribution 
of the occupancy of queue k is governed only by the distribution of and is a 
function of only Afc. Further, we assume that the first packet of every probe pair, 
packet Pi, 'sees' each of the queues in steady state. This allows us to express the 
distribution of d as a function of only A. We then use the results of the previous 
sections to compute the analytical distribution 'I'd (A) for a given value of A. Let 

I>^d,*d(A),^ denote a suitably defined distance between ^ d and 'I'd (A). The 
estimate A of the parameter vector can be obtained as 



A = argmin 



i.e., A is obtained as the minimizer of the distance between the empirical and the 
analytical distributions of the separation vector. 

A well known measure of the distance between two distributions is the Kullback- 
Leibler (KL) distance. The KL distance between the empirical and the analytical 

distributions for a given A, 2?kl ^d(A), *d) , is defined as 

M***-< A >) -E^xiog^^. 

It can be shown that the minimizer of 2?kl (&d> ^d(A)^ is the maximum-likelihood 
estimate of A. Alternatively, we could use the squared Euclidean distance, as the 
measure of the distance between the two distributions. The squared Euclidean 
distance T>e (^d, ^ , d(A), ^ between the empirical and the analytical distributions 
for a given A is given by 

V E (* d) * d (A)) := (* d = x - * d = x ( A )) 2 ' 

X 

A simplistic way to obtain A is as follows. The analytical distribution is com- 
puted and stored for all the values of the parameter vector A over a suitably dis- 
cretized grid on the feasible range of A. We obtain the empirical distribution ^ d 
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from a sufficiently large number of probes. The minimizcr of X>(#d,*d(A)J is 
then obtained by an exhaustive search over the grid. To find a reasonably accurate 
estimate, the exhaustive search will require a very fine grid and consequently a large 
number of distribution calculations. Rather than an exhaustive search, we can also 

use an optimization method to reach the minimum of T> (j&d, 'I'd(A)^ . 

As we had mentioned earlier, for the case of one queue, if we fix do = 1, then 
d\ will have the same distribution as an d estimating the parameters of 
is straightforward. We will not pursue that any more. The case of a two-queue 
path has interesting numerical properties and we explore them in detail in the next 
section. This will give us an insight into the difficulties of probing paths that have 
multiple bottlenecks. We expand upon this in Section [7j We will also see that the 
estimation quality when probing a multicast tree is significantly better. 



6 Numerical Results 

In this section we present some numerical results obtained as follows. Different 
paths and multicast trees of independent discrete-time queues are simulated. The 
probe packets are injected into the system at very low probing rate. The passage 
of the probes through the different queues are simulated and the output separation 
obtained for each probe. We do not simulate lost probes although it is easy to 
see that they do not affect the algorithm. Although our results are applicable for 
all distributions of the arrival process, we report our results for the case when the 
packet arrivals form an i.i.d. Poisson sequence. Since Poisson is a one-parameter 
distribution, it is easier to discuss the issues in the estimation of A. 

We first consider a two-queue path. Let Ai and A2 be the arrival rates of packets 
to queues 1 and 2 respectively. Probe pairs are generated according to a Poisson 
process. We obtain the empirical distribution from 2 x 10 8 slots of simulation time. 
Table Q] shows the estimates that minimize the squared Euclidean distance and 
KL distance between the empirical distribution and the calculated distribution for 
different true values of Ai and A2 and for two probing rates. Notice that both the 
distance measures give comparable precision in almost all cases. Since the squared 
Euclidean distance is computationally simpler than the KL distance, the rest of the 
numerical results are based on the squared Euclidean distance. 

The exhaustive search method has the drawback that it is primarily a one time 
estimation and not an adaptive one. An alternative is to reach the optimum value 
by a suitable optimization algorithm, for example, a steepest descent algorithm. 
This may not however give the best estimate if the distance function has a local 
minima. We performed an extensive numerical investigation of the nature of the 
cost function for the two-queue case. However, the cost function was seen to be 
non-convex. The level sets of the Euclidian cost function as shown in Fig. [6] are 
clearly not convex. This may present some problem in convergence of an optimiza- 
tion technique. However, the limited numerical simulations performed by us gave 
reasonably good performance under ideal setup. 
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Table 1: Exhaustive search based estimates for the cascade of two queues 



True 
values 
(Ai,A 2 ) 


inter 
probe 
mean 


Estimates in 3 different executions 


Euclidean cost based 


KL distance based 


(Ai,A 2 ) 


(Ai,A 2 ) 


(Ai,A 2 ) 


(Ai,A 2 ) 


(Ai,A 2 ) 


(Ai,A 2 ) 


(0.9, 0.8) 


1000 


(0.83,0.85) 


(0.88,0.81) 


(0.87,0.82) 


(0.87,0.82) 


(0.88,0.81) 


(0.88,0.81) 




200 


(0.91,0.79) 


(0.91,0.79) 


(0.87,0.82) 


(0.87,0.82) 


(0.88,0.81) 


(0.88,0.81) 


(0.5, 0.6) 


1000 


(0.54,0.57) 


(0.51,0.59) 


(0.51,0.59) 


(0.54,0.57) 


(0.51,0.59) 


(0.50,0.60) 


JJ 


200 


(0.47,0.62) 


(0.51,0.59) 


(0.49,0.61) 


(0.51,0.59) 


(0.50,0.60) 


(0.49,0.61) 


(0.7, 0.4) 


1000 


(0.59,0.50) 


(0.71,0.39) 


(0.74,0.36) 


(0.68,0.42) 


(0.71,0.40) 


(0.68,0.42) 


55 


200 


(0.70,0.40) 


(0.69,0.41) 


(0.69,0.41) 


(0.70,0.40) 


(0.71,0.39) 


(0.70,0.40) 


(0.2, 0.3) 


1000 


(0.22,0.28) 


(0.13,0.36) 


(0.35,0.15) 


(0.22,0.28) 


(0.13,0.36) 


(0.21,0.29) 


5! 


200 


(0.20,0.30) 


(0.19,0.31) 


(0.21,0.29) 


(0.20,0.30) 


(0.20,0.30) 


(0.20,0.30) 



0.6^ 




1 

\hat{\lambda_2} 

\hat{lambda_1) 



Figure 5: A typical Euclidean cost as a function of (Ai, A 2 ). This plot was 
obtained with Ai = 0.3 and A 2 = 0.8. 
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Figure 6: Contours of the Euclidean cost function between the analytical 
and the empirical distributions for Ai = 0.3 and A2 = 0.8. The contours 
were plotted by taking logarithm of the cost function so that more details 
can be seen near the valley (that almost touches zero). 
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Fig. [5] shows the plot of a typical Euclidean cost function for an empirical 
distribution obtained from simulation with (Ai,A2) = (0.3,0.8). The minima is at 
(0.31,0.79). However, notice that the function shows a distinct valley along a line 
with an extremely low slope along it. The existence of the valley is more clear from 
the contour plot of the cost function shown in Fig. [S] The steepest descent method 
will result in a very slow convergence due to zig-zagging along the valley. 

Our main simulation involves an iterative stochastic gradient based online es- 
timation. We divide the probe sequence into equal sized blocks and compute the 
'instantaneous' empirical distribution of d for every block. In our simulations we 
have used a block size of 500. Let ^Jf^, be the instantaneous empirical distribution 
obtained from the n-the block of probes. The overall empirical distribution after 
the 7j-th iteration, ^d.n, is updated using the exponential update rule 



Here a is a constant between and 1. We used a = 0.05 in all our simulations. 
For each n, we then update our estimate of A based on the current empirical 
distribution ^d,n- We use the steepest descent algorithm for one iteration for each 
n and update the estimate of A as 



where | ■ | denotes the length of a vector and a n is the step size. We used an adaptive 
step size so that whenever a particular step increases the distance, the step size is 
decreased and whenever the gradient is more than a threshold, the step size is 
increased. The increase and decrease of the step size is not done if it is respectively 
above or below some thresholds. If the arrival process is non stationary, this scheme 
allows us to track slowly varying A. Our simulations will study the effectiveness of 
such tracking. 

We first consider the two-queue path. Ai and A2 are varied linearly with time 
and a step change is also simultaneously introduced. Fig. [7] plots the true values 
of Ai and A2 and the estimates as functions of time. For the estimation, the initial 
estimates are taken as Ai = A 2 = 0.5. From these plots we see that the steepest 
descent algorithm converges and tracks the true values reasonably well even when 
there is step change. It is instructive to study the evolution of the estimate. Fig. [S] 
shows the evolution of the estimate pair (Ai, A2). A close look at the figure reveals 
that initially as well as when the parameters change suddenly, the estimates go to 
the nearest point on the new valley and then zig-zags along the valley towards the 
true values. 

We next consider the three-queue multicast tree of Fig. |4] As with the two- 
queue case, the true values of Ai are varied linearly and also a step change is 
introduced. In Fig. [9] we plot Ai, A 2 and A3 as functions of time along with the true 




(13) 



A n = A 
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Figure 7: The estimates of Ai and A2 and their true values for two-queue 
path. The straight lines represent the true values of the parameters. 
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Figure 8: The evolution of the estimates of Ai and A2 and their true values 
for the two-queue path. The straight lines show how the true value of (Ai, A2) 
changes. 

values. Observe that these estimates are far superior to those for the two-queue 
path. This is because here we have the empirical joint distribution of the output 
separations at two leaf nodes of the network and this gives much more information 
about the arrivals to the queues. 

The slow convergence of the adaptive estimation due to the valley in the Eu- 
clidean cost function that we saw for the two-queue path is expected to be even 
more pronounced for a three-queue path. In fact, in our simulations, the estimates 
did not converge to the true values even after 2 x 10 8 slots of simulation time. We 
conjecture that there will be a valley along a two dimensional manifold. 

We have seen that in a multicast network, having multiple outputs through 
different paths helps improve the quality of the estimation. However two or more 
queues in series makes for slow convergence of the estimates. To investigate the 
interplay of these two we carried out the estimation for a four-queue multicast tree 
shown in Fig. llOl The plots of the estimates and the true values of the arrival means 
are shown in Fig. Illl As expected, we observed faster convergence for this network. 
However, the queues 2 and 4 are in series in this network. The estimates of A2 and 
A4 can be seen to be worse than the estimates for the two-queue path. This is 
possibly because of the error in the estimate of Ai (this misleads the estimation of 
A2 and A4 with an incorrect distribution of the probe separation at the input of 
queue 2) and also because the probe pairs entering queue 2 are not separated by a 
fixed one slot. Hence, the convergence of the estimates of A2 and A4 is not as fast 
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Figure 9: The estimates of Ax, A2 and A3 and their true values for the 
multicast network 
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Figure 10: The four-queue multicast tree used in some of the simulation 
results. 
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Figure 11: The estimates of Ai, A2, A3 and A4 and their true values for the 
network of Fig. [TUl 

as that of Ai and A3. 

In the adaptive estimation, we need to consider two performance criteria. If 
the arrivals are stationary, we need to minimize the error in the estimate. However, 
if the arrivals are non-stationary, we would like the estimates to track the true 
values closely. Achieving these objectives is governed by the choices of the various 
parameters of the estimation algorithm. In the rest of this section we will discuss 
the role of these parameters. The performance of the estimation was seen to be 
quite sensitive to the selection of the parameters, and quite a bit of tuning was 
necessary to get the performance shown in the figures. 

Choice of a, block length and a n : The parameter a (in Eqn. [T3|) and the number 
of probe pairs in a block dictate how fast the empirical distribution tracks the 
true distribution of the output separations. However, there is a trade-off between 
tracking the true distribution for non-stationary arrivals and the variance of the 
empirical distribution under stationary arrivals. Increasing a or decreasing the 
block size enables the empirical distribution to track the actual distribution faster 
but increases the variance of the empirical distribution under stationary arrivals. 
This reflects in A„ tracking the true value faster but in the components of A„ 
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having larger variance under stationary arrivals. 

The step size a n also plays a role in the convergence rate and the tracking 
ability of the estimation. A larger a n enables faster tracking but results in larger 
variance of the estimates under stationary arrivals. 

How well to track the minima of V (&d,n, (A)^ : In our simulations, wc have 
used one iteration of steepest descent for each n. Alternatively we could use an 
exhaustive search to obtain the minima of T> (^&d.n, (A)^ . As a second alter- 
native, we may execute many iterations of steepest descent for each n. This may 
involve a fixed number of iterations or termination based on a specified condition. 
A larger number of iterations for each n will enable A„ to be nearer the minima of 

the current cost function D (&d,n, 'J'd(A)^ . Hence, a larger number of iterations 

for each n will result in the estimate A„ tracking the minima of V (*d,n,*d(A)) 
more closely. Although this results in better tracking, it makes the estimate more 
tuned to the current empirical distribution and as a result more sensitive to the 
variation of the current empirical distribution under stationary arrivals. 
Alternative Optimization Methods: The valley in the cost function makes the con- 
vergence of steepest descent algorithm very slow due to the typical zig-zag path 
followed by this algorithm. Other optimization algorithms like conjugate gradient 
or Newton's algorithm may also be tried to overcome the problem. 
Method of Moments Based Estimation: A common method in parameter estimation 
is the method of moments. A similar approach can be considered here. We can 
compute a few moments of the calculated and the empirical distribution and try to 
match them by varying the parameters. In the absence of closed form expressions for 
the moments, we may again need to estimate by minimizing a suitable cost function 
between the calculated moments and the moments from the empirical distributions. 
However, since we do not have a way of calculating the moments directly without 
calculating the distributions first, this method presents little attraction, especially 
since we can expect to lose some information in working only with some moments 
instead of the distributions themselves. Indeed, the results obtained by using the 
squared Euclidean distance between the vectors of the first ten moments of the 
calculated and the empirical distributions gave no promising results. It is also to 
be noted that computation of higher order moments is computationally intensive 
and needs numerical caution in dealing with large numbers. 

Probing rate: Finally, we remark on the probing rate. A higher probing rate gives us 
more samples and hence a quicker convergence for stationary arrivals and a better 
ability to track fast changing arrival processes. However, the probing process makes 
the aggregate arrival process, including the probe packets, non Poisson and disturbs 
the estimate because the calculated distribution will not be correct. Further, the 
samples may not be independent if the probes are closely spaced. To minimize the 
effect, the mean probe arrivals should be small. In our simulations, the times at 
which probes P\ are inserted were generated according to a Poisson process of rate 
0.005. 
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Figure 12: The convention for A n , X n and Q n used in Section [7,11 

7 An Extension and a Generalization 

In this section we consider an extension and a generalization to the algorithms 
presented in Sections [3] and 2J First, we relax the assumption that the probe packet 
has the least priority within its batch. We will see that this essentially involves 
extension of the joint distribution of A\ tTn and X^ m _ x of Section ETT1 to include 
packets that arrive in the same slot as the probes but may be queued after the 
probe packet. We then present a generalization in which we allow the transfer 
delays in the network to be random rather than fixed but arbitrary. 

7.1 Generalization to Equal Priority for the probes 

In our discussions earlier, we assumed that the probe packets have the lowest pri- 
ority among the packets that arrive in the slot. We now consider the case when 
all the packets that arrive in a slot have the same priority and the position of the 
probe packets within the batch may have a distribution. 

For convenience, we redefine the notation and let Q n be the number of packets 
in the queue at the beginning of slot n but before the arrival of the packets in slot n. 
This new notation is shown in Figure [T2l Under this notation the queue evolution 
described by Eqn. Q] changes to 

Qn+l = Qn + A n — X n , 

Qn+m Qn F ^Fi,n+m— 1 -^n.n+m— 1 • (1^) 

Note that according to this definition the stationary distribution of Q n would 
be different from that given in Eqn. [2J The probability generating function for the 
stationary distribution would be 

aw . li^fcii 

z — V\z) 
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and the recursion for the probabilities will be 



l-A 
Po 

£K-i(i-pi)- 

TTQPg) 



for g = 
for g = 1 

for g > 1. 



For a probe packet arriving into the queue in slot k, define A' k (resp. A/-) to 
be the number of packets from those that arrived in slot k that are queued before 
(resp. after) the probe packet. Thus A' k + Ak = Ak- Also, for probe packets arriving 
in slots i and j, define, 



j' 

i,, 



A 



h3 



A, 
A, 
A, 



A l+1 
A l+ i 
A i+ i 



Aj-i 



4 

A', 



As before, we assume that the probe packets Pi and P 2 enter the queue in slots 
and di respectively. They depart in slots D\ and D 2 . Along the lines of Eqn. 0] 
we obtain the output separation of the probe packets as follows. 



D 2 



= Qo + A' Q 



di ■ 
di ■ 
D 2 

di - 



A' 



l + A 



O.di 



Qo 

-Dx 

1 + A a ^i — Xo,di- 



X, 



O.di 



A' 



(15) 



For any di 



1 + 1 > 1, to obtain the distribution of the output separation, we 
thus need to know the joint distribution of Ao >m and ^4o,m+i- 

The distribution of A' m and A m is derived from the distribution of Aq as 

00 

Pr (A' m = k) = Pr (I m =k)=J2 J^yPr (Ao - j) ■ 



Following the procedure in obtaining the recursion in Section[¥]we begin by defining 
Xo m to be the number of departures in slots 0, . . . , m when Qo — q and Aq = a. 
As before, we begin by assuming Qq = Aq = and compute the joint distribu- 



tion of Aq'^ and A[ m . The base equation will be 



Pr X, 



(x° o: =l,A lfi =j) 



1 for I — 1 and j 
otherwise. 



(16) 



Here ^4i,o is an empty summation and we define it to be 0. We first obtain the 
recursion for the joint distribution of Ag'^ and A\, m and then the joint distribution 
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of Xq'^ and A' l m+1 . Following the arguments of Section[3l we obtain the following 
recursions. 



Pr 





0,0 



l,m-l — / — 1 

xPr(A m =0) 6{j-l+l) 

xPr(A m =j-n), 

for Z < to + 1, and j > I — 2, 
otherwise 



(17) 



j 

Pr - I, A 1>m =j-k) Pr (A' m+1 = k) , 



(18) 



k=0 



where S(-) is the Kronecker delta function. If Qq = q and Aq = a, then, following 
the arguments in the development of Eqn. 1101 the joint distribution of X^^ and 
A'i m+1 would be obtained in exactly the same manner as Eqn. 1 1 01 except that q 
would be replaced by q + a. Thus we will have 



P r v^o]m ~ ^i,m+i — 3) 

q — a. A' l m+1 = j ^ 
for I < to + 1 and q + a < I — 1 , 



t= m +l 



Pr X, 



(19) 



+ i „ ■' ' yii.m ~ *>^i,m+l — 3j 
for Z = m + 1 and q + a < / — 1 , 
Pr(A' 1>ro =j) 

for Z = m + 1 and q + a > Z — 1 , 
otherwise. 

Of the a packets that arrive in slot 0, if a are queued after Pi, then we need that 
many fewer arrivals in 1, . . . , to. Hence, we have 



-1=3 



1 



(20) 



Finally, we can uncondition on q and a as 
Pr (Xo.m — I, Ao, m+ i = j 



EE P r(^ 

q=0 a=0 

Pr (Q = g) Pr = /, A , ro+1 = j) 



(21) 
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Figure 13: The multicast tree with random transfer delays. 



The output separation of the probes can now be obtained following the same 
method as in Section [3] 

7.2 Random Delays on the Links 

In deriving the distribution of the output separation, in Eqns. [TT] and 1121 we as- 
sumed that the transfer delays from the output of a queue to the input of the next 
queue on the path is arbitrary but fixed. We can relax this requirement as follows. 

Consider the series of queues of Figure [3] first. In going from the output of 
queue k to the input of queue k + 1, let Pi be delayed by t\ slots and P2 by t% slots. 
Let dk be the probe separation at the output of queue k. The probe separation at 
the input to queue k + 1 is d' k := dk + (t% — It is reasonable to assume that the 
order of the probes is maintained in the transfer from k to k + 1. Hence t\ and t 2 
cannot be independent. If the conditional distribution of d' k given dk is known then 
Eqn. [TT] can be generalized as 

00 

Pr(4-!=i) = £Pr(d*-i=Jkfc-i=i) 

i=l 

00 

Pr (dk = s) = Pr (d k = s|4-l = 3) Pr (4-1 = j) 
J=i 

Along similar lines, we can generalize for random transfer delays on the multi- 
cast tree. The schematic is shown in Fig. [T3l Eqn. [T2l can be generalized as 

Pr (d 2 = r,d 3 = s) 

OO OO OO 

71— 1 m— 1 l—l 

x Pr (tii, 2 = I] tii, 3 = Tn\d\ = n) Pr (d\ = n) 

There are many uses for this generalization. If the transmission rates on the 
links are different, then the separation is compressed or expanded depending on 
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whether the downstream link has a higher or lower data rate respectively. The 
conditional distribution of the separation can be obtained by knowing the two 
transmission rates. 

A second use is in tomography of paths with many links with one or two of them 
being 'bottleneck links.' Consider a if -queue path. Let links i and j, 1 < i < j < K , 
be the bottleneck links. The delays from the source to the input of link i, from the 
output of i to input of j and from the output of j to the destination can have simple 
models. For example, let K = 3, i = 1 and j — 3 with link 2 being lightly loaded. 
We can assume that over link 2, each of the probes Pi and P2 will see at most one 
cross traffic packet independently. Let a be the probability of this event. Then we 
can see that, if the input separation is greater than 1, the separation increases by 
1 if only P2 is sees a cross traffic packet, while it shrinks by 1 if only Pi sees such 
a packet. If the input separation is 1, then the output separation will also be 1 if 
neither P2 does not see a cross traffic packet. We can then use 



Pr (g?2 = n\di = m) 

a if 771 — 1 and n = 2 

1 — a if m = n = 1 

a(l — a) if n = 77i+l or n = m — 1 & m > 1 

1 — 2a(l — a) if ri = m and m > 1 

otherwise. 



Here link 2 is playing the role of the random delay element. 



8 Discussion and Conclusion 

We conclude with a discussion on the connection with some earlier results on the 
analysis of packet-dispersion based probing and also on network delay tomography. 

We first explore the connection of this work with that in [211 US] ■ A theoretical 
analysis of packet pair probing of a single bottleneck path is carried out in [21] . 
A continuous time model is used and the key result is Theorem 2. This theorem 
states that for a probe pair inserted into the queue at time ai with input separation 
S and probe packet length s, the output separation 5' is given by 

*' = ^ + ^ + max(M^£ ) (22) 

where 

• C is the capacity of the bottleneck link, 

• Ys(i) is the work (sum of the service times of the packets) that enters the 
queue in (t, t + S) and 

• Bg(t) is the unused capacity of the server in the interval (t, t + S). 



29 



Thus the output separation is a sample of a linear combination of three interdepen- 
dent random processes. This theorem essentially implies that obtaining the joint 
law of these processes is the key to the use of packet pair probing. 

We can interpret Eqn. [2H for the discrete time queues that we consider as 
follows. Note that s = C = 1. The first term is the number of packets entering the 
queue between the enqueuing of the two probe packets and the third term is the 
number of idle slots between the enqueuing of the two probe packets. Thus we can 
rewrite Eqn. [22] as follows. 

dout = Ai idi +l + max{<i i -Xo,d i -i,0} 
= Ai t di + 1 + di — Xo,d 4 _i. 

The last equality is true because the number of departures in di slots can never be 
greater than di. This is nothing but Eqn. 0) 

Eqn. [22] essentially says that the output separation of a packet pair is a function 
of the work entering the queue between the probe arrivals and the amount of work 
done in this interval. These two quantities arc dependent and one needs their joint 
distribution for finding the distribution of d ou t . The key difficulty in developing the 
analytical model for packet pair probing is to obtain this joint distribution, and we 
believe, that is a key contribution of this paper. 

It is interesting to contrast our problem setting with delay tomography of, for 
example, 8, 9], since both estimate how the individual links will affect the traffic 
flowing, albeit, in different forms. An important assumption in the development of 
network-delay tomography is that the end-to-end delays of the individual packets 
can be measured. This requires that the sources and the destinations be suit- 
ably time synchronized. This is clearly not an easy task. In fact, the packet-pair 
techniques were developed to avoid this problem. Furthermore, much of the delay- 
tomography results were developed assuming that the probes are multicast packets. 
This means that when there are links in series, the estimation is for the delay across 
the series and not on the individual links. 

Finally, we remark that the primary objective in this paper is to develop an 
analytical framework for 'capacity tomogrpahy' and not develop a practical tool. 
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